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ABSTRACT 


Highly  maneuverable  ocean  going  vehicles  require  quick  response,  control,  and 
guidance  to  maintain  accurate  track  keeping  characteristics.  Ocean  vehicles,  however,  may 
experience  significant  lags  in  their  inertial  positional  information  that  limit  their  reaction 
times.  This  thesis  investigates  the  effects  of  these  lags  on  guidance  and  control.  The 
relationship  of  critical  visibility  versus  the  controller  time  constant  and  its  effect  on  the 
stability  of  the  guidance/control  scheme  is  analyzed.  Results  are  presented  based  on  time 
domain  and  frequency  response  techniques  using  a  dynamic  model  of  the  Naval 
Postgraduate  School  Autonomous  Underwater  Vehicle  II  (NPS  AUV  II),  for  which  a 
complete  set  of  hydrodynamic  coefficients  and  geometric  properties  is  available. 
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I .  INTRODUCTION 


A.  BACKGROUND 

Marine  vehicles  that  require  high  maneuverability,  such  as 
hydrofoils,  require  quick  response,  control  and  guidance  to 
maintain  accurate  track  keeping  characteristics.  This  is 
accomplished  through  successful  path  planning,  navigation, 
guidance  and  autopilot  design  [Ref.  1]. 

Sufficient  information  is  obtained  from  charted  obstacles 
and  the  environment  for  smooth  paths  to  be  generated  for  the 
vehicle  to  follow  [Ref.  2].  A  certain  level  of  feedback  is 
provided  through  the  use  of  sonar  and  acoustics  in  order  to 
replan  a  path  when  uncharted  obstacles  are  encountered  or  when 
the  mission  objectives  are  changed.  Based  on  inertial 
positional  information,  the  guidance  law  provides  the 
appropriate  vehicle  heading  by  suitable  use  of  the  vehicle 
effectors  such  as  rudders,  dive  planes,  and  cross  body 
thrusters.  However,  lags  in  obtaining  and  processing  inertial 
positional  information  can  limit  vehicle  reaction  time.  In 
addition,  the  guidance  law  must  be  as  fast  as  possible  in 
order  to  ensure  accurate  path  keeping  characteristics  [Ref. 
3]  .  Therefore,  the  stability  of  the  combined  guidance/control 
scheme  becomes  an  issue  that  needs  to  be  analyzed. 
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B.  OBJECTIVE  OF  THIS  THESIS 


This  thesis  investigates  the  effects  of  positional 
information  time  lags  on  guidance  and  control.  The 
relationship  of  the  critical  visibility  (minimum  vehicle 
lookahead  distance)  versus  the  controller  time  constant  and 
its  effect  on  the  stability  of  the  combined  guidance/control 
scheme  is  analyzed.  Results  are  presented  based  on 
computations  using  time  domain  and  frequency  response 
techniques.  All  computations  are  performed  for  a  dynamic 
model  of  the  NPS  AUV  II  [Refs.  4  and  5]  for  which  a  complete 
set  of  hydrodynamic  coefficients  and  geometric  properties  is 
available . 

C.  THESIS  OUTLINE 

Chapter  II  presents  the  vehicle  equations  of  motion  in  the 
horizontal  plane,  the  equations  that  govern  steering  control, 
and  the  guidance  scheme  used  in  this  research. 

Chapter  III  presents  the  stability  analysis  based  on 
eigenvalue  computation  of  the  combined  guidance  and  control 
scheme  developed  in  Chapter  II.  Results  are  presented  based 
on  the  stability  curves  developed  for  first,  second  and  third 
order  approximations  for  time  lag  in  the  commanded  vehicle 
heading  in  the  control  law. 

Chapter  IV  presents  an  exact  computation  of  the  stability 
curves  based  on  frequency  response  methods,  which  utilizes  the 
Nyquist  criterion  for  stability. 
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II.  EQUATIONS  OF  MOTION,  CONTROLLER  DESIGN  AND  GUIDANCE  SCHEME 

A.  INTRODUCTION 

The  equations  of  motion  that  describe  maneuvering  of  a 
marine  vehicle  in  the  horizontal  plane  are  presented  in  this 
chapter.  A  linear  full  state  steering  feedback  control  law  is 
designed  based  on  the  three  equations  in  sway,  yaw,  and  rate 
of  change  of  heading  angle.  The  guidance  scheme  is  based  on 
pure  pursuit  guidance  for  path  following  along  straight  line 
segments . 

B.  EQUATIONS  OF  MOTION 

Restricting  our  attention  to  motion  in  the  horizontal 
plane  (steering  control),  the  mathematical  model  consists  of 
the  nonlinear  sway  and  yaw  equations  of  motion  only.  With  a 
moving  coordinate  frame  fixed  at  the  vehicle's  geometric 
center,  Newton's  equations  of  motion  are 

miv+ur+x^-YoT2) -Y  (2.1) 

Izi:+mxG(v+ur) -myGvr=N  (2.2) 
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where  v  and  r  are  the  relative  sway  and  yaw  velocities  of  the 
moving  vehicle  with  respect  to  the  water;  m  is  the  mass  of  the 
vehicle;  y^,  yG  are  the  respective  lateral  and  longitudinal 
positions  of  the  center  of  gravity;  and  Y,N  represent  the 
total  excitation  sway  force  and  yaw  moment,  respectively. 
These  forces  can  be  expressed  as  the  sum  of  quadratic  drag 
terms  and  first  order  memoryless  polynomials  in  v  and  r,  which 
represent  the  added  mass  and  damping  due  to  the  vehicle's 
motion  through  the  water.  In  this  way,  Y  and  N  can  be 
represented  by 


Y  =  -|l4y,,/+-|t3(y^+rriir)  +  -|  t.2Yvuv 


nose 

tail 


5) 


Iv'+Srl  2  0 


N  =  sNti+£p  (N*v+Nrur)  +£(iNvuv 
2  2  2 


where  l  is  the  vehicle  length,  and  Ya,  Na  represent  partial 
derivatives  of  Y  and  N  with  respect  to  a,  is  the  drag 

coefficient,  and  8  is  the  rudder  angle. 

Equations  (2.1)  and  (2.2)  can  be  nondimensionalized  with 
respect  to  the  constant  forward  speed  u,  and  the  vehicle 
length  {  ;  with  the  dimensionless  time  variable  being  tu/l  . 
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The  cross  flow  integral  drag  terms  become  important  for 
hoveling  operations  or  low  speed  manuevering,  whereas  at 
relatively  high  speeds  and  low  angles  of  attack  (with  respect 
to  the  water),  the  steering  response  is  predominantly  linear. 

The  model  becomes  complete  with  the  addition  of  the 
expressions  for  the  vehicle  yaw  rate 


(2.3) 


and  the  inertial  position  rate  (horizontal  plane) 


y=usin«Jr+vcost|f  (2.4) 

where  is  the  heading  angle  as  shown  in  Figure  1. 

Equations  (2.1),  (2.2),  and  (2.3)  can  be  written  as  a  set 

of  three  coupled  linear  differential  equations  of  the  form 


ijr=r 


(2.5a) 


v=a^uv+a12iir+biU2b 


(2.5b) 


i =a12uv+a22ur+b2u2& 


(2.5c) 
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where  the  coefficients  and  hx  are  functions  of  the  vehicle 
hydrodynamic  and  geometric  properties,  y  is  the  heading  angle 
and  6  is  the  rudder  angle. 

The  linearized  form  of  equation  (2.4)  at  a  nominal  u  is 


y=LnJr  +  v  (2.6) 

Hence,  the  final  set  of  linearized  equations  of  motion  used 
for  steering  control  are 


(2.7a) 


v^a11uv+a12ur+Jb1u2fi  (2.7b) 


r =a12uv+a22ur+Jb2u2fi  (2.7c) 


y=ui|/+v  (2 . 7d) 


at  any  nominal  u. 

This  system  of  equations  becomes  the  basis  for  the 
controller  design. 
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C .  CONTROLLER  DESIGN 

Equations  (2.7a),  (2.7b),  and  (2.7c)  govern  the  steering 

control  of  the  model  used  in  this  research. 

These  equations  can  be  represented  in  the  form 


x=Ax+bt> 


(2.8) 


where  the  state  vector  equation  is 

x  =  [ijr,  v,  r]  T  (2.9) 

Linear  full  state  feedback  is  introduced  in  the  form 

6-k^  ( i|r-tjr0)  +k2v+k3r  (2.10) 

where  \\fc  is  the  commanded  heading  angle,  and  the  gains  kt,  k2, 
k3  are  computed  by  pole  placement  such  that  the  closed  loop 
system  of  equations  (2.8),  (2.10)  has  the  desired  dynamics. 

The  existence  of  the  difference  (\|/-yc)  in  the  control  law 
(2.9)  has  the  effect  of  trying  to  point  the  vehicle's 
longitudinal  axis  towards  the  desired  heading. 

If  the  desired  characteristic  equation  is  selected  so  that 
its  time  constant  is  tc  dimensionless  seconds,  its  general 
form  becomes 
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(2.11) 


A,3+a1A.2+«2A.+a3  =  0 


where 


{bxa22-b2ax2)  uJk2  +  {b2alx-b^a2X)  u%k^a2+b2u2k^  (2.12b) 

Jb1u2k2+i32u2ic3=-a1-  (a11+a22)  u  (2  •  12c) 

Due  to  the  explicit  dependence  on  u,  the  gains  kj,  k2,  k3  are 
continuously  updated  for  different  nominal  forward  speeds. 

D .  GUIDANCE  SCHEME 

The  guidance  scheme  used  in  this  research  is  an 
orientation  based  command  scheme.  In  this  scheme,  the 
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commanded  vehicle  heading  angle  yc  is  a  function  of  the  line 
of  sight  angle  a  between  the  actual  vehicle  position  and  a 
reference  position  on  the  nominal  path  which  is  located  at  a 
constant  lookahead  or  visibility  distance  d  ahead  of  the 
vehicle  as  shown  in  Figure  1. 

This  line  of  sight  angle  is  defined  by 


tano  = 


y 

d 


(2.13) 


The  autopilot  is  then  called  upon  to  deliver  the  commanded 
heading  angle  \|JC .  The  simplest  orientation  guidance  law  is 
pure  pursuit  guidance  where  the  commanded  heading  angle  \(fc  in 
the  control  law  (2.9)  equals  the  line  of  sight  angle  a  in 
(2.13)  . 

Then  yc  is  defined  by 


t|rc  =  -arc  tan  ^ 


(2.14) 


For  relatively  small  angles 


=  -arctan-^  » 
c  d  d 


(2.15) 
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This  results  in  a  three  state  heading  autopilot  design 
which  can  exhibit  very  robust  characteristics.  However,  the 
stability  of  this  scheme  may  become  an  issue  when  transiting 
along  straight  line  segments  since  the  commanded  heading  angle 
is  a  function  of  the  vehicle  response,  i.e.,  autopilot 
response  is  limited  by  time  lags  in  vehicle  dynamics.  We  will 
next  examine  the  stability  of  this  scheme. 
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III.  EIGENVALUE  ANALYSIS 


A.  INTRODUCTION 

In  this  chapter,  we  present  the  stability  analysis,  based 
on  eigenvalue  computation,  of  the  combined  guidance  and 
control  scheme  developed  in  Chapter  II.  A  brief  presentation 
of  the  stability  curve  is  developed  first  for  the  case  of  zero 
time  lag.  Incorporation  of  a  nonzero  time  lag  in  the 
positional  information  formally  results  in  an  infinite  state 
space  system.  This  is  truncated  into  first,  second,  and  third 
order  approximations.  The  resulting  eigenvalue  problems  are 
then  analyzed  with  emphasis  on  stability  curves. 


B.  STABILITY  CONSIDERATIONS  AND  TIME  LAG 

1.  Stability  Considerations 

In  pure  pursuit  guidance,  where  the  commanded  heading 
angle  \|fc  in  the  control  law  equals  the  line  of  sight  angle  o, 
the  trivial  equilibrium  state  which  corresponds  to  straight 
line  motion  is  characterized  by 


i|f  =  v=r=y=0 


(3.1) 


Linearization  of  the  state  equations  in  the  vicinity  of  (3.1) 
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produces  the  linear  system 


x-Ax, 


(3.2) 


where  the  complete  state  vector  is 


x-  [ilr ,  v,  r,y]  T 


(3.3) 


Local  stability  properties  of  (3.1)  are  then  established  by 
the  eigenvalues  of  A.  The  characteristic  equation  of  (3.2)  is 
a  quartic  of  the  form 


X4  +BX3+CX2+DX+E=0 , 


(3.4) 


where  the  coefficients  B,  C,  D,  E  are  functions  of  the  vehicle 
properties,  the  lookahead  distance  d,  and  the  controller  gains 
k1#  k2,  k3.  A  pair  of  complex  conjugate  roots  of  (3.4)  crosses 
the  imaginary  axis  when  the  following  conditions  are  met 
[Ref .  6] . 


bcd-b2e-d2=o, 


(3.5a) 


(3.5b) 
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a  (^Ia22  -^2ai2  ^2^  ^lg1't~  (^ia22~.^2ai2~-^2^  C3  jj  Yc) 

{b2axl-b1a21)2u 

The  conditions  (3.5a)  and  (3.5b)  determine  the  critical 
visibility  dcric  for  stability.  For  d  >  d^,.,  the  equilibrium 
state  (3.1)  is  stable  which  means  that  the  control  law  will 
drive  and  keep  the  vehicle  onto  the  straight  line  path.  For 
d  <  dcirit,  the  equilibrium  state  loses  its  stability  and  the 
vehicle  response  becomes  oscillatory  as  a  result  of  the  pair 
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of  complex  conjugate  eigenvalues  with  positive  real  parts. 

Results  of  the  critical  visibility  versus  the 
controller  time  constant  t_  for  zero  time  lag  are  presented  in 
Figure  2.  As  expected,  the  higher  values  of  tc  (i.e., softer 
controller)  require  higher  lookahead  distances  d  for  path 
accuracy.  This  agrees  with  existing  guidelines  that  the 
guidance  law  must  be  sufficiently  slower  than  the  controller 
in  order  for  the  dynamics  of  the  two  not  to  interfere  with 
each  other  and,  therefore,  guarantee  stability  of  the  combined 
scheme.  Very  high  values  of  d  correspond  to  a  very  slow 
guidance  law  with  a  loss  in  speed  of  response  and  path 
accuracy  [Ref.  1]. 

2 .  Time  Lag 

In  ocean  vehicles,  all  states  needed  in  the  control 
law  are  readily  available  at  the  required  rate,  with  the 
possible  exception  of  the  positional  information  y  (Figure  1)  . 
The  latter  may  require  a  significant  data  analysis  and 
reduction  of  sonar  returns  and  inertial  navigational 
information.  As  a  result,  it  is  likely  that  a  time  lag  will 
exist  in  the  positional  information  y  which  is  used  in  the 
guidance  law. 

With  the  introduction  of  a  time  lag  T  (sec),  the  commanded 
rudder  angle  \|fc  in  the  pursuit  guidance  control  law  becomes 
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0  0.5  1  1.5  2  2.5  3 


d 

Figure  2.  Critical  Visibility  Distance  d  versus  te  for 
Zero  Time  Lag  T. 
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*„  =  -aictani^l 
c  d 


(3.8) 


and  for  relatively  small  angles 


y(t-T) 

d 


(3.9) 


The  resultant  rudder  control  law  (2.9)  becomes 


fi  =k^  +k1  +k2  v+k3r 


(3.10) 


After  some  algebra,  the  resultant  linearized  equations  of 
motion  (2.7a),  (2.7b),  (2.7c),  (2.7d)  become 


i|r=r  (3.11a) 


<'»2>1u2Je1ilr+  (a11u+i?1u2ic2)  v+fa^u+^u2^)  r+i71uz-^y(  t-T)  ( 3  .  lib) 


f-=ix!u2£1i|i+  ( a21u+b2u2k2 )  v+  ( a2Zu+b2u2k3 )  r+b2u2-~y{  t-T)  (3 . 11c) 
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y=UTjr  +  V 


(3 .lid) 


The  Taylor  expansion  for  the  term  y(t-T)  is 


yU-T)  « 


T’Z  'T'i 

y~Ty+~-y-—y+.  .  . 


(3.12) 


We  can  now  examine  the  stability  of  the  previously  developed 
control  scheme  for  first,  second  and  third  order 
approximations  for  time  lag  in  the  control  law. 


C.  FIRST  ORDER  APPROXIMATION  FOR  TIME  LAG 

For  a  first  order  approximation  for  time  lag  T 


y(  t-T)  m  y-Ty 


(3.13) 


where 


$ - ui|r + v'  (3.14) 

After  some  algebra,  the  resultant  linearized  equations  of 
motion  (3.11a),  (3.11b),  (3.11c),  (3. lid)  become 
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(3.15a) 


t=r 


y=  (b1u2Jc1-l)1u3— 1 T)  t|r+  (a11u+Jb1u2k2-.b1u: 


D  v+ (a^u+bju2^)  r+i1u2-^y(  3 


15b) 


ic  ic  ic 

f =  (i32u2^1-Jb2u3-^r)  T|r+  {a21u*b2u2k2-b2u2~  T)  v+  ( a22u*b2u2k2 )  i+b2u2-^y^  3  ^5C  j 


y=injr  +  v 


(3 .15d) 


The  state  space  form  of  the  equations  of  motion  is 


x-Ax 


(3.16) 


In  matrix  form,  (3.15a),  (3.15b), 


(3.15c),  (3.15d)  become 


♦ 

0 

0 

l 

0 

V 

(a3lu*b1u2k2-b1u2-^T) 

(aijU+^u2*,) 

t 

(£2u2*l-£au3-^n 

(a21u+b2u2k3-b2u2-^T) 

(a„u+.ft2uJ*j) 

y 

u 

l 

0 

0 

* 


V 

r 


(3.17) 


y 


The  standard  eigenvalue  problem  (3.17)  is  solved  to 
determine  numerically  the  critical  d  versus  tc  curve  for  a 
given  value  of  T.  Appendix  A  presents  the  program  used  for 
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the  first  order  approximation  for  time  lag. 

Figure  3  shows  the  resulting  stability  curves  for  the 
first  order  case.  It  shows  the  relation  between  the  critical 
visibility  d  and  the  controller  time  constant  tc  for  values  of 
time  lag  of  0.0,  0.5,  1.0,  1.5  and  2.0  seconds. 

As  expected,  the  stability  curve  shifts  to  the  right  with 
increasing  values  of  time  lag.  For  the  same  time  constant  tc, 
the  critical  visibility  distance  increases'  with  corresponding 
increases  in  the  amount  of  time  lag. 

D.  SECOND  ORDER  APPROXIMATION  FOR  TIMS  LAG 

For  a  second  order  approximation  for  time  lag  T,  we  have 


T2 

y(t-T)  =  y-Ty+±-? 


(3.18) 


where 


y=in|f+v,  y=ur+ v  (3.14),  (3.19) 


Following  some  algebra,  the  resultant  linearized  equations  of 
motion  (3.11a),  (3.11b),  (3.11c),  (3. lid)  become 
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Figure  3.  First  Order  Approximation  For  Time  Lag: 

Critical  Visibility  Distance  d  versus  tc  for 
Time  Lags  T  =  0.0,  0.5,  1.0,  1.5,  2.0  sec . 
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(3.20a) 


where 


♦  -r 


v=A2iilr  +A22v+A23r+A2iy 


(3.20b) 


r=A3it  +Aj  2  v+A3  3  r  +A34y 


(3.20c) 


y=ui|r  +  v 


(3 .20d) 


•^21  = 


Jb1u2ic1-jb1u3-^T 

i-2> lUdr2 

1  2d 


■^22  = 


an  u  +bz  u  zk2  — Jbx  u  2  ~  r 

1-i,  U2^r2 

1  2d 


•^23  = 


a12u+b1u2jc3+b1u3-|^:r2 

1-i,  U2_^l  r2 

1  2d 
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(b2u2^) +(b2u2^±T2) 


(b^u2-!) 

(l-D^Ar2) 

20 


As  with  the  first  order  case,  the  state  space  equations  of 
motion  reduce  to 
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x-Ax 


(3.16) 


In  matrix  form,  (3.20a),  (3.20b),  (3.20c),  (3.20d) 


* 

0  0  10 

* 

V 

■^21  "^22  -^23  ^24 

V 

i 

■^31  A32  -^33  ^34 

T 

y 

U  1  0  0 

y 

become 


(3.21) 


Simular  to  the  first  order  case,  the  standard  eigenvalue 
problem  (3.21)  is  solved  using  Fortran  programming  and  MATLAB 
in  order  to  develop  the  stability  characteristics  of  the 
second  order  case.  Appendix  B  presents  the  program  used  to 
develop  the  stability  characteristics  for  this  case. 

Figure  4  shows  the  resulting  stability  curves  for  this 
case.  The  results  of  the  critical  visibility  d  versus  the 
controller  time  constant  tc  for  time  lags  of  0.0,  0.5,  1.0, 
1.5  and  2.0  secs  are  shown. 

As  with  the  first  order  case,  the  higher  values  of  tc 
require  higher  lookahead  distances  d  for  path  accuracy.  High 
values  of  d  correspond  to  a  slow  guidance  law  with  a  loss  in 
speed  of  response  and  path  accuracy,  and  the  stability  curve 
shifts  to  the  right  with  increasing  values  of  time  lag. 
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Figure  4.  Second  Order  Approximation  For  Time  Lag: 

Critical  Visibility  Distance  d  versus  tc  for 
Time  Lags  T  =  0.0,  0.5,  1.0,  1.5,  2.0  sec. 
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However,  due  to  being  more  accurate  than  the  first  order  case, 
the  second  order  approximation  provides  tighter  stability 
curves  for  the  same  time  lags. 

E.  THIRD  ORDER  APPROXIMATION  FOR  TIME  LAG 

For  a  third  order  approximation  for  time  lag  T 


y(t-T) 


2  <t*3 

y-Tt*—9-—y 


(3.22) 


where 


y=ui|r  +  v,  y-uz+v,  y-  uz + v  (3.14),  (3.19),  (3.23) 

After  a  considerable  amount  of  algebra  and  following  a  simular 
procedure  as  with  both  the  previous  two  cases,  the  resultant 
linearized  equations  of  motion  reduce  to  the  state  space  form 

Bx-Ax  (3.24) 


or 


x=B~1Ax 
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The 
(3.11b) , 


where 


resultant  linearized  equations  of  motion 
(3.11c),  (3. lid)  become 


t  *r 


B„r+B35v2  =  A^+A22vsA^v2+A2ir+A^ 


B53i+B55v2  -  -A51Tjr  +A52 vx  +A53 v2  +A54r  +AS5y 


y=ur^+v2 


v,  =  v. 


B 


33 


A 
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2?1u2/Ci-iJ1u3 


d 


T 


(3.11a) , 


(3.25a) 


(3.25b) 


(3.25c) 


( 3 . 25d) 


( 3 . 25e) 
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A32  =  a11u+£1u2/c2-Jb1u2-±:r 


£ 

A33  =  a12u+b3u2k3+bxuz 


35 


and 


B 


53 


=  l+jb2u3-^T3 
2  6  c? 


sS5  =  *»*£»- 


*31 


b2u2-b2u3 


h 

d 


T 


ki, 


as2  =  a21u+b2uzk2-b2u2  -±T 
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,3  *1 


^53  =  a22u+i32u2ic3+i)2u3^r2 


A54  =  b*u2^j 
ASS  =  b2u2^jT2 

In  matrix  form,  (3.25a),  (3.25b),  (3.25c),  (3.25d),  (3.25e) 

become 


1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

* 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

*33 

0 

*35 

t 

- 

Aj  1 

*32 

*33 

*3* 

*35 

r 

0 

0 

0 

1 

0 

y 

u 

1 

0 

0 

0 

y 

0 

0 

*53 

0 

*55 

A>i 

*52 

*53 

*54 

*55 

Simular  to  both  the  first  order  and  second  order  cases, 
the  generalized  eigenvalue  problem  (3.21)  is  solved  using 
Fortran  programming  and  MATLAB  in  order  to  develop  the 
stability  characteristics  for  the  third  order  case.  Appendix 
C  presents  the  program  used  to  develop  the  stability 
characteristics  for  this  case. 

Figure  5  shows  the  resulting  stability  curves  for  this 
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d 

Figure  5.  Third  Order  Approximation  For  Time  Lag: 

Critical  Visibility  Distance  d  versus  te  for 
Time  Lags  T  =  0.0,  0.5,  1.0,  1.5,  2.0  sec. 
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case.  The  results  of  the  critical  visibility  d  versus  the 
controller  time  constant  tc  for  time  lags  of  0.0,  0.5,  1.0, 
1.5  and  2.0  seconds  are  shown.  As  with  the  first  and  second 
order  cases,  the  higher  values  of  tc  require  higher  lookahead 
distances  d  for  path  accuracy,  high  values  of  d  correspond  to 
a  slow  guidance  law  with  a  loss  in  speed  of  response  and  path 
accuracy,  and  the  stability  curve  shifts  to  the  right  with 
increasing  values  of  time  lag.  However,  due  to  being  more 
accurate  than  either  of  the  other  two  cases,  the  third  order 
approximation  provides  the  tightest  stability  curves  for  the 
same  time  lags. 


F.  COMPARISON  OF  RESULTS 

Results  of  the  critical  visibility  d  versus  the  controller 
time  constant  tc  for  first,  second  and  third  order 
approximations  for  time  lag  have  been  obtained.  Figures  3,  4 
and  5  presented  the  stability  curves  for  first,  second  and 
third  order  approximations  for  time  lags  of  0.0,  0.5,  1.0,  1.5 
and  2.0  secs.  These  figures  have  shown  that  stability  of  the 
control  scheme  decreases  with  increasing  time  lag,  and  for  the 
same  time  constant,  the  critical  visibility  distance  increases 
with  increasing  time  lag. 

Figures  6,  7,  8  and  9  present  a  comparison  of  first, 
second  and  third  order  approximations  for  each  time  lag 
considered.  These  figures  show  that  for  each  time  lag,  the 
third  order  model  presents  the  largest  region  of  stability  for 


31 


Critical  Visibility  Distance  d  versus  tc  for 
First,  Second  and  Third  Order  Approximation 
Cases  and  Time  Lag  T  =  0.5  sec. 


Figure  7.  Critical  Visibility  Distance  d  versus  te  for 
First,  Second  and  Third  Order  Approximation 
Cases  and  Time  Lag  T  =  1.0  sec . 
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Critical  Visibility  Distance  d  versus  tc  for 
First,  Second  and  Third  Order  Approximation 
Cases  and  Time  Lag  T  =  1.5  sec . 


Critical  Visibility  Distance  d  versus  tc  for 
First,  Second  and  Third  Order  Approximation 
Cases  and  Time  Lag  T  =  2.0  sec . 


vehicle  straight  line  path  accuracy.  In  addition,  for  the 
same  time  constant,  the  critical  visibility  is  least  for  the 
third  order  case.  This  demonstrates  that  the  minimum 
lookahead  distance  required  for  straight  line  motion  stability 
decreases  with  controller  time  lag  accuracy. 

It  can  be  seen  that  the  first  order  approximation  is  the 
most  conservative  for  design  purposes  since  it  predicts  the 
highest  required  minimum  visibility  distance.  As  expected, 
the  differences  among  the  three  approximations  are  more 
pronounced  as  the  time  lag  increases  and  as  the  controller 
time  constant  becomes  smaller;  i.e.,  tighter  control. 

Figures  10,  11,  12  and  13  present  the  differences  in 
critical  visibilities  among  first,  second  and  third  order 
approximations  for  each  individual  time  lag  considered.  A 
comparison  of  these  figures  demonstrates  that  the  difference 
in  critical  visibility  between  respective  orders  (time  lag 
models)  decreases  with  decreasing  time  lag.  This  indicates 
that  for  low  time  lags,  controller  time  lag  accuracy  is  not  as 
crucial  for  maintaining  vehicle  straight  line  path  accuracy  as 
that  of  higher  time  lags.  The  greater  the  time  lag,  the 
greater  the  difference  in  critical  visibilities  among  time  lag 
models,  and  the  more  crucial  is  the  controller  time  lag 
accuracy  for  stability. 

In  general,  if  the  guidance  law  is  slow  enough  to  allow 
sufficient  time  for  the  controller  to  react,  stability  of 
straight  line  motion  is  guaranteed  [Ref.  1]. 
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dT-dO 


Figure  10.  Difference  of  Critical  Visibility  Distances  at 
Time  Lags  of  T  =  0 . 5  and  0.0  sec  versus  tc  for 
First,  Second  and  Third  Order  Cases. 


Figure  11.  Difference  of  Critical  Visibility  Distances  at 
Time  Lags  of  T  =  1.0  and  0.0  sec  versus  te  for 
First,  Second  and  Third  Order  Cases. 


o  O 


.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 

dT-dO 

Figure  12.  Difference  of  Critical  Visibility  Distances  at 


Time  Lags  of  T  =  1.5  and  0 . 0  sec  versus  te  for 
First,  Second  and  Third  Order  Cases. 


39 


oO 


dT-dO 

Figure  13.  Difference  of  Critical  Visibility  Distances  at 
Time  Lags  of  T  =  2.0  and  0 . 0  sec  versus  te  for 
First,  Second  and  Third  Order  Cases. 
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6.  NORMALIZATION  OF  RESULTS 

In  an  attempt  to  characterize  the  stability  analysis 
results  for  all  values  of  time  lag,  we  proceed  as  follows:  If 
the  vehicle  side  slip  velocity  v  is  very  small  so  that  it  can 
be  neglected,  the  equations  of  motion  (2.7a),  (2.7b),  (2.7c), 
and  (2.7d),  take  the  form 


i|r=r 


t=a22ur+b2u2b 


y=ui|r 


since  v  =  0.  The  control  law  is 

+k*r' 

where  the  first  order  approximation  for  the  commanded  heading 
angle  is 


y-Ty 

~~d~’ 


Based  on  the  above  equations,  we  can  get  the  characteristic 
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equation  of  the  system  as 


s3- {a22u+b2k2u2)  s2- (b2kxuz-b2kru*T^ )  s-bjlqu3-^  =  0. 


Applying  Routh's  criterion  to  this  cubic  equation,  we  find 
that  loss  of  stability  occurs  when 


( a22u+b2k3u 2)  2 -b^k^u 3 T^)  =  -b2kru2^, 

from  which  we  can  get  the  critical  visibility  distance 


d  =  - - —  +  Tu 

a22+bzk2u 


Equation  (3.27)  can  be  written  as 


dr'dc  . 
T 


U, 


(3.27) 


(3.28) 


where  d0  is  the  critical  visibility  distance  for  T  =  0,  and  dj. 
the  critical  visibility  distance  for  a  nonzero  T.  In 
dimensionless  quantities,  equation  (3.28)  is 
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(3.29) 


^r~^o  _  2 
T 


where  u  =  2.0  ft/sec  is  the  vehicle  speed,  and  C  =  7.3  ft  is 
the  vehicle  length. 

A  comparison  among  the  approximate  expression  (3.29)  and 
the  first,  second,  and  third  order  approximations  is  shown  in 
figures  14,  15  and  16.  The  agreement  is  better  for  the  first 
order  approximation,  as  expected,  and  also  for  the  higher 
controller  time  constant  tc.  This  is  because  a  high  tc  results 
in  soft  vehicle  response  with  negligible  amounts  of  side  slip. 
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Figure  16.  Normalized  Critical  visibility  versus  te  for 
Third  Order  Approximation  For  Time  Lags 
T  »  0.5,  1.0,  1.5,  2.0  secs. 
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IV.  FREQUENCY  RESPONSE  ANALYSIS 


A.  INTRODUCTION 

The  previous  analysis  using  Taylor  expansions  of  y(t-T) 
breaks  down  for  approximations  beyond  third  order.  The  reason 
for  this  is  that  for  higher  than  third  order  expansions,  the 
matrix  B  in  the  generalized  eigenvalue  problem  (3.23)  becomes 
singular.  Therefore,  it  was  felt  that  a  different  technique 
should  be  sought  in  order  to  obtain  an  exact  computation  of 
the  stability  curves  and  also  to  check  the  validity  of  our 
calculations.  In  this  chapter,  we  present  a  technique  based 
on  frequency  response  methods,  which  utilizes  the  Nyquist 
criterion  for  stability.  The  resulting  equation  is  solved 
using  Newton's  iteration. 

B.  NYQUIST  CRITERION 

From  our  previously  developed  model,  the  linearized 
horizontal  equations  of  motion  are 


i '=r 


(2.7a) 


v=a11uv+a12u.r+271u2d  (2.7b) 
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±=a12uv+a22ur+b2u2& 


(2.7c; 


y=u\ |r  +  v 


( 2 . 7d) 


at  any  nominal  u.  The  control  law  is 


8=&x  (t-i|»c)  +k2v+k3r 


(2.10) 


where  the  commanded  heading  angle  is 


—  y(fc-D 


(3.9) 


The  Laplace  transform  of  equation  (3.9)  is 


i|rc  =  -4e"Ts 


(4.1) 


and  the  resultant  rudder  control  law  (2.10)  becomes 


lr 

_ —  x/a  ~Ts 


8  = jqijf  +k2v+k3r+-±ye 


(4.2) 


Following  some  algebra,  the  resultant  linearized  equations  of 
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motion  (2.7a),  (2.7b),  (2.7c),  [2.16.)  become 


t=r 


v^(b1uJic1-2)1u3-^Dil»+(a11u+i31u2ic2-2)1uJ-^D  v+  (a^u+i^u2*,)  r+hjU^ye 


ir  k  k 

t=  U>2u2ic1-i>2u3-^n  ♦  +  {.a21u+b2u2k2-b2u2  -~T)  v+  [a22u+b2u2k2)  r+^^-^ye 


y=u\ |r  +  v 


Equations  (4.3a),  (4.3b),  (4.3c),  (4.3d)  reduce  to 

space  form 


x=Ax 


The  matrix  form  becomes: 


♦ 

0 

0 

1 

0 

y 

(i^u2^) 

(a^u+jt^u2^) 

[aX2u*b2u2k%) 

[b1u2~e'TB) 

r 

[b2u3k2) 

(a21u+£2u2/c2) 

(a22u+2J2u2&3) 

(b2u2^e-T') 

y 

u 

1 

0 

0 

(4.3a) 


®>  (4.3b) 


r.  (4.3c) 


(4.3d) 


the  state 


(3.16) 


♦ 

(4.4) 

r 

y 
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The  characteristic  equation  of  the  form  [A-Is]  =  0  becomes 
after  a  considerable  amount  of  algebra 


s4+A3s3+A2s2+A1s+  [B2s2+B1s+B0)  De~Ta  =  0  (4.5) 


where 


A3  =  -  {a11-*-a22)  u  -  {bxk2+b2kz)  u2 


A2  =  (^33^22  ^12^21^  U2+  (^11^2  ^21^1 )  ^ 3 -^3  +  ^22^1  ^12^2^  ^k2~b2U2k.i 


Aj  =  -  (,a12b2  a22bx)  u  k-jDa  Ts-  (a21Jb1~ai1b2)  u^k^ 


B2  =  -i^u2^ 


B1  = 


=  ^b2-a2xbx)  u4iq 


and 


Z?  = 


1 

d 
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The  characteristic  equation  of  the  system  is  written  as 


1  +  K  G  ( s)  =  0 


(4.6) 


where 


(B^+BjS+Bt)  e~Ts 
s  ( s  3  +  A3  s  2  +  Aj  s +  ) 


(4.7) 


is  the  transfer  function,  and  we  have  denoted 

K  -  D  (4.8) 


With  s=j(0,  the  phase  angle  <j>  is  represented  by 


$  =  ^G(JU) 


(4.9) 


and  it  follows  that 


<|>  =  Z(e'7yw)  -  Z(j'w)  +  L  (-B2ci)2+fl1ju+fl0)  -  L  (-jw3-A3Ci>2+A2jca+Al) 


4>  =  -wr 


n 

2 


+ 


tan'1  ( 


Bi<*  ) 
B0-B2u>2 


tan'1  ( 


-u3+A,&)  . 

- -£_) 

-A3(i)2+A1 


(4.10) 
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For  stability,  the  Nyquist  criterion  states  (Figure  17)  that 
at  the  solution  of  the  equation 


4>(0>)  =  -71 


(4.11) 


which  is  phase  cross  over  frequency  CO  =  G^,  the  gain  margin 
must  be  equal  to  1 


\KGljtoi)  |  =  1 


(4.12) 


where  K  =  D  =  1/d.  It  follows  that  equation  (4.11)  becomes 


+tan~1  ( — — -1  ■ 

2 


)  -tan'1  ( 


1)  =  —%  (4.13) 


Rearranging,  we  get 


tan*1  ( 


B0-B2(o 


i 


-tan"1  ( 


-AjCOi+A,^ 


-n  +  — +w,T 
2  1 


Taking  the  tangent  of  both  sides  and  using  the  identity 


tan(x-y) 


_  tan  {x) -tan  (y) 
1+tan  (x)  tan  (y) 
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Three  Different  Values  of  Gain  [Ref.  7]. 
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and  rearranging  following  some  algebra,  the  result  becomes 


-(50-Ba(P i)  _  l 

(B0-B2Qi)  (-Ajtal+Aj  +S1«1(-w?+A2&)1)  tan(«1T) 


From  equation  (4.12),  it  follows  that 


1  \-Baal+B1jul+Bii _ 


Solving  for  d,  we  get 


J  (B0-B2<*l)2+B?ul 

(^-AjOi)  2+  (AjWj-Wi)  2 


With  a  time  lag  T  =  0,  we  get  from  equation  (4.14) 


(B0-B2«i)  +  B1(01  (-Ui+AjQj)  =0 


Rearranging,  we  get 


(ByAj-BjJtt?  ♦  (B1A2-B2A1-B0A3)  o»i  +  B^  =  0 


(4.14) 


(4.15) 


(4.16) 
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Equation  (4.16)  can  be  solved  exactly  for  o  and  then  d 
can  be  computed  from  (4.15).  This  zero  time  lag  calculation 
was  performed  in  order  to  validate  the  results  of  d  versus  tc 
for  the  zero  time  lag  solution  of  equation  (3.6)  .  The  results 
using  the  two  different  techniques  were  found  to  be  identical. 

C .  ALGORITHM  DESCRIPTION 

Introducing  the  terms 

oc  i  -  BZA 3  ~B-± 

fit  2  =  B1A2  -B2&L  ~BqAj 

*3  =  *0*1 

and 

Pi  *  -*2 

P  3  =  *1*1  ~*0*2 

P2  *  *0 -*-*2*2 -*1*3 
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Then  equation  (4.12)  becomes 


P1W5+P20)3+P3£«)  _  1 

a1(i>4+a2ci>2+a3  tan(<oD 


Rearranging,  we  get 


(P1&>5+P2ci)3+P3ci>)  tan  (tor)  +  a1o4+a2ci)2  +  a3  » 


Equation  (4.18)  is  now  in  the  form 


f  (to)  =  0 


Following  Newton's  iteration  for  0),  we  have 


=  <■>*-! 


*'<•*- 1> 


From  equation  (4.20), 

■f(u)  =  (p1ci)5+P2Q3+P3ti))  sin(wD 
+  (a1Q4+a2u2+a3)  cos  («r) 


(4.17) 


(4.18) 


(4.19) 


(4.20) 


(4.21) 
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and 


f'iui)  =  (5P1o>4+3P2ci>2+P3)  sin(G>r) 

+  (P1G)5+P2<1>2  +  P,<I>)  TCOS(UiT) 

{A. 22) 

+  (4a1o>3+2a2o)  cos  (uT) 

-  (a1G)4+«2o2+a3)  rsin(o)T’) 

Fortran  programming  and  MATLAB  were  used  in  order  to 
develop  the  stability  characteristics  for  the  Newton's 
iteration  method.  Appendix  D  presents  the  program  used  to 
develop  the  stability  characteristics  for  this  case. 

Figure  18  shows  the  resulting  stability  curves  for  the 
Newton's  iteration  method.  The  results  of  the  critical 
visibility  d  versus  the  controller  time  constant  tc  for  time 
lags  of  0.0,  0.5,  1.0,  1.5,  and  2.0  secs  are  shown. 

As  with  the  first,  second  and  third  order  approximation 
cases,  the  higher  values  of  tc  require  higher  lookahead 
distances  d  for  path  accuracy,  high  values  of  d  correspond  to 
a  slow  guidance  law  with  loss  in  speed  of  response  and  path 
accuracy,  and  the  stability  curves  shift  to  the  right  with 
increasing  values  of  time  lag.  However,  since  the  Newton's 
iteration  method  presents  the  exact  solution  of  the  frequency 
response  to  equation  (4.19),  the  resulting  stability  curves 
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d 

Figure  18.  Newton's  Iteration  Method:  Critical  Visibility 
Distance  d  versus  te  for  Time  Lags  T  =  0.0, 
0.5,  1.0,  2.0  secs. 
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represent  the  precise  stability  characteristics  for  these  time 
lags . 


D.  RESULTS  AND  DISCUSSION 

Results  of  the  critical  visibility  d  versus  the  controller 
time  constant  tc  for  the  Newton's  iteration  method  has  been 
demonstrated.  Figure  18  presented  the  precise  stability 
curves  for  time  lags  of  0.0,  0.5,  1.0,  1.5  and  2.0  secs. 

This  figure,  as  well  as  those  for  the  first,  second  and 
third  order  approximation  cases,  demonstrates  that  the 
stability  of  the  control  scheme  decreases  with  increase  in 
time  lag  and  for  the  same  time  constant,  critical  visibility 
increases  with  increase  in  time  lag. 

Figures  19,  20,  21,  22  and  23  present  a  comparison  of 
generated  stability  curves  among  the  Newton's  iteration  method 
and  those  for  the  first,  second  and  third  order  approximation 
cases  for  time  lags  of  0.0,  0.5,  1.0,  1.5  and  2.0  secs.  It 
can  be  seen  that  there  is  barely  a  slight  difference  in  the 
stability  curves  for  the  Newton's  iteration  method  case  and 
the  third  order  approximation  case.  This  indicates  that  a 
third  order  approximation  for  time  lag  presents  the  best  model 
for  accounting  for  a  positional  information  time  lag  in  the 
vehicle  control  law. 
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0  0.5  1  1.5  2  2.5  3 

d 

Figure  19.  Critical  Visibility  Distance  d  versus  tc  for 
Newton's  Iteration  Method  and  First,  Second  and 
Third  Order  Approximations;  T  =  0.0  sec. 
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d 

Figure  20.  Critical  Visibility  Distance  d  versus  te  for 
Newton's  Iteration  Method  and  First,  Second  and 
Third  Order  Approximations;  T  =  0.5  sec. 
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d 

Figure  21.  Critical  Visibility  Distance  d  versus  tc  for 
Newton's  Iteration  Method  and  First,  Second  and 
Third  Order  Approximations;  T  =  1.0  sec. 
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d 

Figure  22.  Critical  Visibility  Distance  d  versus  tc  for 
Newton's  Iteration  Method  and  First,  Second  and 
Third  Order  Approximations;  T  =  1.5  sec. 
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d  vs  tc  for  1st,  2nd,  3rd  Ord  tc  newtons  method  time  lag  T=2.0 


CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

A  methodology  for  considering  positional  information  time 
lags  in  the  control  law  for  vehicle  maneuvering  in  the 
horizontal  plane  has  been  presented.  The  relationship  of  the 
critical  visibility  versus  the  controller  time  constant  for 
both  zero  and  non-zero  time  lags  has  been  established. 

Time  lags  have  been  approximated  by  first,  second  and 
third  order  Taylor  expansions  in  the  commanded  vehicle  heading 
in  the  control  law.  A  comparison  of  the  resulting  stability 
curves  has  demonstrated  that  the  third  order  approximation 
presents  the  greatest  stability  and  the  least  critical 
visibility  for  a  given  time  lag.  It  was  found  that  the  first 
order  approximation  was  the  most  conservative  for  controller 
design  purposes  since  it  predicted  the  highest  required 
critical  visibility  distance.  In  addition,  the  differences 
among  the  three  approximations  became  more  pronounced  with 
increasing  time  lag  and  decreasing  controller  time  constant 
(tighter  control). 

Further  analysis  was  conducted  to  characterize  the 
stability  results  for  any  value  of  time  lag  using  a 
normalization  procedure  and  the  first  order  approximation  for 
time  lag  in  the  commanded  vehicle  heading  in  the  control  law. 
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A  comparison  of  the  normalized  stability  curves  for  each  time 
lag  showed  that  the  greatest  differences  occur  with  increasing 
order  and  decreasing  controller  time  constant.  This  is  due  to 
a  tight  (quick)  vehicle  response  with  significant  side  slip  at 
lower  time  constants. 

Frequency  response  techniques  based  on  the  Nyquist 
stability  criterion  were  used  to  produce  the  exact  stability 
curves  for  comparison  with  the  previously  generated  Taylor 
expansion  stability  curves.  These  results  virtually  matched 
those  obtained  for  the  third  order  approximation  time  lag  case 
and  validated  our  calculations  and  stability  curves  generated 
for  the  case  of  zero  time  lag.  These  results  indicated  that 
the  third  order  approximation  best  represents  a  positional 
information  time  lag  in  the  control  law,  however  for  design 
purposes,  the  first  order  approximation  provides  the  most 
conservative  approach. 


B .  RECOMMENDATIONS 

Some  recommendations  for  further  research  are  as  follows: 

*  Experimental  verification  using  the  NPS  AUV  II. 

*  Consideration  of  time  lag  in  the  simulation  of  an  Inertial 
Navigation  System  required  for  positional  updates. 

*  Analysis  of  positional  information  time  lags  in  the  motion 
stability  in  the  vertical  plane,  along  curved  paths,  or 
with  respect  to  other  guidance  schemes. 
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APPENDIX  A 


PROGRAM  THESIS1 . FOR  (Time  Delay-lst  Order  Approx) 
BIFURCATION  ANALYSIS 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1 , K2 , K3 , L , NR , NV, NDRS , NDRB , I Z , MASS , 
&  NRDOT, NVDOT 

DIMENSION  A(4,4),FVl(4),IVl(4),ZZZ(4,4),WR(4),WI(4) 

OPEN  ( 11 , FILE-' BI FI .RES' , STATUS- 'NEW ) 

OPEN  (12, FILE- ' BI F2 . RES ' , STATUS- ' NEW’ ) 

OPEN  (13,FILE-'BIF3.RES' , STATUS- ' NEW' ) 

Vehicle  Parameters: 

WEIGHT-435.0 
IZ  -45.0 
L  -7.3 
RHO  -1.94 
G  -32.2 
XG  -0.0104 
MASS  -WEIGHT/G 

YRDOT  — 0.00178*0. 5*RHO*L**4 
YVDOT  — 0.03430*0. 5*RHO*L**3 
YR  -+0. 01107*0. 5*RHO*L**3 
YV  — 0.03096*0. 5*RHO*L**2 
YDRS  -+0. 02345*0. 5*RHO*L**2 
YDRB  -+0. 02345*0. 5*RHO*L**2 
NRDOT  —0.00047*0. 5*RHO*L**5 
NVDOT  — 0. 00178*0. 5*RHO*L**4 
NR  — 0.01022*0. 5*RHO*L**4 
NV  — 0.00769*0. 5*RHO*L**3 
NDRS  — 0. 337*0 .02 345*0. 5*RHO*L** 3 
NDRB  -+0. 283*0.02345*0.5*RHO*L**3 

WRITE  (*,1001) 

READ  (*,*)  TCMIN , TCMAX, ITC 

WRITE  (*,1002) 

READ  (*,*)  XDM I N , XDMAX , I XD 

XDMIN-XDMIN*L 
XDMAX-XDMAX  *  L 
WRITE  (*,1003) 

READ  ( * , *  )  U 

WRITE( * , 1100 ) 

READ  ( * , * )  TL 

DH  - ( I Z -NRDOT) * ( MASS- YVDOT ) - 
&  ( MASS *XG- YRDOT ) * ( MASS  *XG-NVDOT ) 

AAl 1- ( ( IZ-NRDOT ) * YV- ( MAS S*XG- YRDOT ) *NV )/DH 
AA12«(  ( IZ-NRDOT)  M-MASS  +  YR)- 
&  ( MASS *XG- YRDOT ) * ( -MASS  *XG+NR ) )/DH 

AA2 1- ( ( MASS- YVDOT ) *NV- ( MASS  *XG-NVDOT ) *  YV ) /DH 
AA22-( ( MASS -YVDOT ) * ( -MASS*XG+NR) - 
&  ( MASS*XG-NVDOT ) * ( -MASS+YR ) )/DH 

BB 1 1 - ( ( I Z-NRDOT ) *  YDRS- ( MASS  *XG- YRDOT ) *NDRS ) /DH 
BBl2-( ( IZ-NRDOT ) *YDRB- ( MAS S*XG- YRDOT) *NDRB )/DH 
BB21-( (MASS- YVDOT) *NDRS-(MASS*XG-NVDOT)» YDRS )/DH 
BB22- ( (MASS -YVDOT ) *NDRB- < MASS*XG-NVDOT ) * YDRB )/DH 
C 

WRITE  (*,1004) 
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READ  ( * , * ) 


RATIO 


BBl-BBll+RATIO*BBl2 

BB2-BB21+RATIO*BB22 

EPS  -l.D-5 
ILMAX-1500 

DO  1  I-l.ITC 

WRITE  (*,2001)  I , ITC 

TC-TCMIN+ ( 1-1 ) * ( TCMAX-TCMIN )/( ITC-1 ) 

TCD-TC*L/0 
POLEC-1 . 0/TCD 
ADl-3 . 0* POL EC 
AD2-3 . 0*POLEC**2 
AD3-POLEC*  *  3 
Al-BBl *U*U 
Bl-BB2*U*U 

Cl— ADl-(AAll+AA22)*U 
A2- ( BBl *AA2  2-BB2  *AAl 2 ) *U*  *  3 
B2-( BB2  *AAl 1-BB1 *AA21 ) *U*  *  3 
Kl-AD3/( ( BB2*AAl 1-BBl *AA21 ) *U**3 ) 

C2-AD2- ( AAl 1 ‘AA22-AA1 2  *AA2 1 )*U**2  +  BB2*U*U*Kl 
K2-(C1*B2-C2*B1 )/(Al*B2-A2*Bl ) 

K3-(C2*A1-C1*A2)/(A1*B2-A2*B1 ) 

DO  2  J-1,IXD 

XD-XDMIN+ ( J-l ) *(XDMAX-XDMIN)/( IXD-1 ) 

CALL  LINEAR.(  Kl ,  K2  ,  K  3  ,  U ,  XD ,  AA1 1 ,  AAl  2  ,  AA21 ,  AA2  2  ,  BBl ,  BB2  ,  A,  TL ) 
CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL ( DEOS , WR , WI , FREQ ) 

IF  (J.GT.l)  GO  TO  10 

DEOSOO-DEOS 

XDOO  -XD 

LL-0 

GO  TO  2 

DEOSNN-DEOS 

XDNN  -XD 

PR-DEOSNN*DEOSOO 

IF  (PR.GT.0.D0)  GO  TO  3 

LL-LL+1 

IF  (LL.GT.3)  STOP  1000 
I L-0 

XDO-XDOO 

XDN-XDNN 

DEOSO-DEOSOO 

DEO SN- DEOS NN 

XDL-XDO 

XDR-XDN 

DEOSL-DEOSO 

DEOSR-DEOSN 

XD- ( XDL+XDR )/2 . DO 

CALL  LINEAR) Kl , K2 , K3 , U, XD , AAl 1 , AAl 2 , AA21 , AA22 , BBl , BB2 , A, TL ) 
CALL  RG( 4 , 4 , A, WR ,WI,0,ZZZ,IV1, FV1 , I ERR ) 

CALL  DSTABL ( DEOS , WR, WI , FREQ ) 

DEOSN-DEOS 

XDM-XD 


PRL-DEOSL*DEOSM 

PRR-DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

XDO-XDL 

XDN-XDM 

DEOSO-DEOSL 

DEOSN-DEOSM 

IL-IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 
DI F-DABS ( XDL-XDM ) 

IF  (DIF. GT . EPS )  GO  TO  6 

XD-XDM 

GO  TO  4 

5  IF  (PRR.GT.0.D0)  STOP  3200 

XDO-XDM 
XDN-XDR 
DEOSO-DEOSM 
DEOSN-DEOSR 
IL-IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 
DI F-DABS ( XDM-XDR ) 

IF  (DIF. GT . EPS )  GO  TO  6 
XD-XDM 

4  LLL-10+LL 

WRITE  ( LLL , * )  XD/L.TC 
3  XDOO-XDNN 

DEOSOO-DEOSNN 
2  CONTINUE 
1  CONTINUE 


1001  FORMAT 

1002  FORMAT 

1003  FORMAT 

1004  FORMAT 
1100  FORMAT 
2001  FORMAT 

END 


('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF 
('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF 
( '  ENTER  U' ) 

('  ENTER  BOW/STERN  RUDDER  RATIO') 

(  '  ENTER  TIME  LAG  TL' ) 

(215) 


C 


SUBROUTINE  DSTABL( DEOS , WR ,WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR( 4 ) , WI ( 4 ) 

DEOS--1 . 0D+20 
DO  1  1-1,4 

IF  (WR( I ) .LT.DEOS)  GO  TO  1 
DEOS-WR( I ) 

IJ-I 

1  CONTINUE 
OMEGA-WI ( I J ) 

OMEGA-DABS ( OMEGA ) 

RETURN 

END 


TC'  ) 
XD'  ) 


SUBROUTINE  LINEAR ( Kl , K2 , K3 , U , XD , AAl 1 , AA12 , AA21 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  K1,K2,K3 
DIMENSION  A( 4 , 4 ) 

A( 1 , 1 ) -0 . 0D0 
A ( 1 , 2 )-0 .0D0 
A( 1 , 3 )-l . 0D0 
A( 1 , 4 )-0 . 0D0 


,  AA22 , BBl , BB2 , A, TL ) 
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A(2,l)-  BBl *U*U*K1*( XD-U*TL )/XD 

A(2,2)-U*(BB1*U*(K2*XD-K1*TL)+AA11*XD)/XD 

A(2, 3)-AAl2*U+BBl*U*U*K3 

A( 2 , 4  )-  BBl *U*U*K1/XD 

A(3,l)-  BB2*U*U*K1*(XD-U*TL)/XD 

A( 3 / 2 ) «U* (BB2*U*( K2  *XD-Kl *TL ) +AA2 1 *XD ) /XD 

A( 3, 3)-AA22*U+BB2*U*U*K3 

A( 3 , 4  ) -  BB2  *U*U*Kl/XD 

A(4,l)-U 

A(  4 , 2 )-l . ODO 

A<4, 3)-0.0D0 

A<4,4)-O.ODO 

RETURN 

END 
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APPENDIX  B 


PROGRAM  THESIS2 . FOR  (Tine  Delay-2nd  Order  Approx) 
BIFURCATION  ANALYSIS 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  Kl , K2 , K3 , L, NR , NV , NDRS , NDRB , I Z , MASS , 
&  NRDOT, NVDOT 

DIMENSION  A( 4 , 4 ) , FVl ( 4 ) , I Vl ( 4 ) , ZZZ(4,4) ,WR(4) ,WI(4) 

OPEN  (11,FILE-'BIF1.RES' , STATUS- ' NEW ' ) 

OPEN  (12,FILE-'BIF2.RES' , STATUS- ' NEW' ) 

OPEN  (13,FILE-'BIF3.RES' , STATUS- ' NEW' ) 


C 

C 

C 


Vehicle  Parameters: 

WEIGHT-435.0 
IZ  -45.0 
L  -7.3 
RHO  -1.94 
G  -32.2 
XG  -0.0104 
MASS  -WEIGHT/G 

YRDOT  — 0.00170*0. 5*RHO*L**4 
YVDOT  --0. 03430*0. 5*RHO*L**3 
YR  -+0. 01187*0. 5*RHO*L**3 
YV  — 0.03896*0. 5*RHO*L**2 
YDRS  -+0. 02345*0. 5*RHO*L**2 
YDRB  -+0. 02345*0. 5*RHO*L**2 
NRDOT  — 0.00047*0. 5*RHO*L**5 
NVDOT  — 0.00178*0. 5*RHO*L**4  . 

NR  — 0.01022*0. 5*RHO*L**4 
NV  — 0.00769*0. 5*RHO*L**3 
ndrs  — 0.337*0.02345*0.5*RHO*L**3 
NDRB  -+0.283*0.02345*0. 5*RHO*L**3 

WRITE  (*,1001) 

READ  (*,*)  TCMIN.TCMAX, ITC 

WRITE  (*,1002) 

READ  (*,*)  XDMIN , XDMAX , IXD 

XDMIN-XDMIN*L 

XDMAX-XDMAX  *  L 

WRITE  (*,1003) 

READ  ( * , *  )  U 

WRITE) * , 1100 ) 

READ  ( * , * )  TL 

DH  -( IZ-NRDOT)*(MASS-YVDOT)- 
&  ( MASS *XG- YRDOT ) * ( MASS  *XG-NVDOT ) 

AAl l-( ( I Z -NRDOT ) * YV- ( MASS  *XG- YRDOT) *NV ) /DH 
AA1 2- ( ( I Z -NRDOT ) * ( -MASS+ YR ) - 
&  ( MASS*XG-YRDOT) * ( -MASS*XG+NR) )/DH 

AA21-( (MASS-YVDOT) *NV-(MASS*XG-NVDOT)*YV)/DH 
AA2  2-( (MASS- YVDOT) * ( -MASS*XG+NR  )  - 
&  ( MASS*XG- NVDOT ) * ( -MASS+ YR ) )/DH 

BBl 1- ( ( IZ-NRDOT) *YDRS- ( MASS*XG- YRDOT) *NDRS)/DH 
BB 1 2- ( ( I Z -NRDOT ) *  YDRB- ( MAS  S  *  XG- YRDOT ) *  NDRB ) /DH 
BB2 1- ( ( MASS-YVDOT ) *NDRS- ( MASS*XG-NVDOT ) * YDRS )/DH 
BB22-( (MASS-YVDOT) *NDRB- ( MASS *XG-NVDOT) * YDRB )/DH 

WRITE  (*,1004) 
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READ  ( * , * ) 


RATIO 


BBl -BBl 1+ RATIO* BBl 2 
BB2-BB2 1 +RATIO*BB22 

EPS  -l.D-5 
ILMAX-1500 

DO  1  I-1,ITC 

WRITE  (*,2001)  I , ITC 

TC-TCMIN+ ( 1-1 ) *(TCMAX-TCNIN)/< ITC-1 ) 

TCD-TC*L/U 

POLEC-1 . 0/TCD 

ADl-3.0*POLEC 

AD2-3 .0*POLEC**2 

AD3-POLEC**3 

Al-BBl*U*U 

Bl-BB2*U*U 

Cl—  AD1-(AA11+AA22)  *U 
A2-(BB1*AA22-BB2*AA12) *U**3 
B2-(BB2*AA11-BB1*AA21 )*U**3 
Kl-AD3/( ( BB2  *AAl 1 -BBl * AA2 1 ) *U*  *  3 ) 

C2-AD2-(AA11*AA22-AA12*AA21 )*U**2+BB2*U*U*K1 
K2-(C1*B2-C2*B1 )/( A1 *B2-A2*Bl ) 

K3-(C2*Al-Cl*A2 )/( A1 *B2-A2  *Bl ) 

DO  2  J-1,IXD 

XD-XDMIN+ ( J— 1 ) * ( XDMAX-XDMIN )/( IXD-1 ) 

CALL  LI NEAR ( Kl , K2 , K3 , U , XD, AA1 1 , AAl 2 , AA21 , AA22 , BBl , BB2 , A, TL ) 
CALL  RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL  DSTABL ( DEOS , WR , HI , FREQ ) 

IF  (J.GT.l)  GO  TO  10 

DEOSOO-DEOS 

XDOO  -XD 

LL-0 

GO  TO  2 

DEOSNN-DEOS 

XDNN  -XD 

PR-DEOSNN*DEOSOO 

IF  (PR.GT.0.D0)  GO  TO  3 

LL-LL+1 

IF  (LL.GT.3)  STOP  1000 
IL-0 

XDO-XDOO 

XDN-XDNN 

DEOSO-DEOSOO 

DEOSN-DEOSNN 

XDL-XDO 

XDR-XDN 

DEOSL-DEOSO 

DEOSR-DEOSN 

XD- ( XDL+XDR ) /2 . DO 

CALL  LI NEAR ( K1,K2,K3,U,XD , AAl 1 , AAl 2 , AA21 , AA22 , BBl , BB2 , A, TL ) 
CALL  RG( 4 , 4 , A, WR ,WI,0,ZZZ,IV1, FV1 , 1  ERR ) 

CALL  DSTABL ( DEOS, WR,WI, FREQ) 

DEOSM-DEOS 

XDM-XD 


PRL-DEOSL*  DEOSM 

PRR-DEOSR*DEOSM 

IF  (PRL.GT.O.DO)  GO  TO  5 

XDO-XDL 

XDN-XDM 

DEOSO-DEOSL 

DEOSN-DEOSM 

IL-IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 
DIF-DABS{ XDL-XDM) 

IF  (DIF.GT.EPS)  GO  TO  6 

XD-XDM 

GO  TO  4 

5  IF  (PRR.GT.0.D0)  STOP  3200 

XDO-XDM 
XDN-XDR 
DEOSO-DEOSM 
DEOSN-DEOSR 
IL-IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 
DI F-DABS ( XDM-XDR ) 

IF  (DIF.GT.EPS)  GO  TO  6 
XD-XDM 

4  LLL-10+LL 

WRITE  ( LLL , * )  XD/L , TC 
3  XDOO-XDNN 

DEOSOO-DEOSNN 
2  CONTINUE 
1  CONTINUE 


1001 

FORMAT 

(  ' 

ENTER 

MIN, 

MAX, 

AND 

INCREMENTS 

OF 

TC'  ) 

1002 

FORMAT 

(  ' 

ENTER 

MIN, 

MAX, 

AND 

INCREMENTS 

or 

XD'  ) 

1003 

FORMAT 

(  ' 

ENTER 

U'  ) 

1004 

FORMAT 

(' 

ENTER 

BOW/STERN 

RUDDER  RATIO' ) 

1100 

FORMAT 

(  ' 

ENTER 

TIME 

LAG 

TL'  ) 

2001 

FORMAT 

(215) 

END 


SUBROUTINE  DSTABL( DEOS,WR,WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 
DIMENSION  WR( 4 ) , WI ( 4 ) 

DEOS--1 . 0D+20 
DO  1  1-1,4 

IF  (WR( I ) .LT.DEOS)  GO  TO  1 
DEOS-WR ( I ) 

IJ-I 

1  CONTINUE 
OMEGA-WK  IJ) 

OMEGA-DABS ( OMEGA ) 

RETURN 

END 


C 

SUBROUTINE  LI NEAR ( Kl , K2 , K3 , U , XD , AAl 1 , AAl 2 , AA21 , AA2 2 , BBl , BB2 , A, TL ) 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl,K2,K3 
DIMENSION  A ( 4 , 4 ) 

A( 1 , 1 ) -0 . 0D0 
A( 1 , 2 ) -0 . 0D0 
A( 1 , 3 ) -1 . 0D0 
A( 1 , 4 ) -0 . 0D0 
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A(2, 1)-<BB1*U*U*K1-BB1*U*U*U*K1*TL/XD)/ 

£  (l-BBl*U*U*Kl*TL*TL/(2.D0*XD) ) 

A( 2 , 2 ) - ( AAl l*U+BBl*U*U*K2-BBl*U*U*Kl *TL/XD )/ 

£  ( l-BBl*U*U*Kl*TL*TL/< 2.D0*XD) ) 

A( 2, 3 )-<AAl2*U+BBl*U*U*K3+BBl*U*U*U*Kl*TL*TL/(2.DO*XD) )/ 
£  ( l-BBl*U*U*Kl*TL*TL/( 2.D0*XD) ) 

A(  2 , 4 ) - ( BB1 *U*U* Kl/XD )/ 

£  (1-BB1*U*U*K1*TL*TL/(2.D0*XD) ) 

A(  3 , 1  )-(BB2*U*U*Kl-BB2*U*U*U*Kl*TL/XD)  + 

£  (BB2*U*U*Kl*TL*TL/(2.D0*XD) )* 

£  ( BBl *U*U*Kl-BBl *U*U*U*  K1 *TL/XD )/ 

£  (1-BB1*U*U*K1*TL*TL/(2.D0*XD) ) 

A( 3,2)-(AA21*U+BB2*U*U*K2-BB2*U*U*Kl*TL/XD)+ 

£  ( BB2  *U*U*Kl *TL*TL/( 2 . DO  *XD ) ) * 

£  ( AA11 *U  +  BBl *U*U*K2-BBl *U*U*  Kl *TL/XD ) / 

£  ( l-BBl*U*U*Kl*TL*TL/( 2 .DO*XD) ) 

A( 3 , 3 ) •( AA22*U+BB2  *U*U*  K3+BB2  *U*U*U*Kl *TL*TL/( 2 . DO  *XD ) )  + 
£  (BB2*U*U*K1*TL*TL/(2.D0*XD) )* 

£  ( AAl 2  *U+BB1 *U*U*K3  +  BB1 *U*U*U*  Kl *TL*TL/( 2 . DO  *XD ) )/ 

£  ( 1-BB1 *U*U*Kl *TL*TL/( 2 . DO  *XD ) ) 

A ( 3, 4 )-( BB2*U*U*Kl/XD)+(BB2*U*U*Kl*TL*TL/( 2 . DO*XD) ) * 

£  ( BBl *U*U*Kl/XD )/ 

£  ( 1-BBl *U*U*Kl *TL*TL/( 2 . DO  *XD ) ) 

A<  4, 1 )-U 
A( 4 , 2 )-l . ODO 
A( 4 , 3 ) -0 . ODO 
A( 4 , 4 ) “0 . ODO 
RETURN 
END 
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no  n  non 


APPENDIX  C 


PROGRAM  THESIS3 . POR  (Time  Delay-3rd  Order  Approx) 
BIFURCATION  ANALYSIS 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DOUBLE  PRECISION  Kl , K2 , K3 , L, NR , NV, NDRS , NDRB , IZ , MASS , 

&  NRDOT , NVDOT 

DIMENSION  A( 5 , 5 ) . B( 5 , 5 ) , FV1 ( 5 ) , I Vl ( 5 ) , ZZZ( 5 , 5 ) , ALFR( 5 ) , 
&  ALFI (5), BETA (5), WR ( 5 ) , WI ( 5 ) 

OPEN  ( 11 , FILE-'BIFl .RES' , STATUS- ' NEW ) 

OPEN  < 12, FILE-'BIF2.RES' , STATUS- 'NEW' ) 

OPEN  (13, FILE- 'BIF3. RES' , STATUS- ' NEW ' ) 

Vehicle  Parameters: 

WEIGHT-435.0 
IZ  -45.0 
L  -7.3 
RHO  -1.94 
G  -32.2 
XG  -0.0104 
MASS  -WEIGHT/G 

YRDOT  --0. 00178*0. 5*RHO*L**4 
YVDOT  — 0.03430*0. 5*RHO*L**3 
YR  -+0. 01187*0. 5*RHO*L**3 
YV  — 0.03896*0. 5*RHO*L**2 
YDRS  -+0. 02345*0. 5*RHO*L**2 
YDRB  -+0. 02345*0. 5*RHO*L**2 
NRDOT  — 0.00047*0. 5*RHO*L**5 
NVDOT  — 0.00178*0. 5*RHO*L**4 
NR  — 0.01022*0. 5*RHO*L**4 
NV  — 0.00769*0. S*RHO*L**3 
NDRS  — 0.337*0. 02345*0. 5*RHO*L**3 
NDRB  -  +  0.283*0.02345*0. 5*RHO*L**  3 

WRITE  (*,1001) 

READ  (*,*)  TCMIN , TCMAX , ITC 

WRITE  (*,1002) 

READ  (*,*)  XDMIN, XDMAX , IXD 

XDMIN-XDMIN*L 

XDMAX-XDMAX*  L 

WRITE  (*,1003) 

READ  (*,*)  U 

WRITE( *,1100) 

READ  ( * , *  )  TL 

DH  -( I Z-NRDOT)*( MASS- YVDOT )- 
6  ( MASS *XG- YRDOT ) * ( MASS  *XG-NVDOT ) 

AAl l-( ( I Z -NRDOT ) * YV- ( MASS  *XG- YRDOT ) *NV )/DH 
AAl 2- j ( I Z -NRDOT ) * ( -MASS+YR ) - 
i  ( MASS *XG- YRDOT ) * ( -MASS*XG+NR ) )/DH 

AA2 1- ( ( MASS- YVDOT ) *NV- ( MASS  *XG-NVDOT ) *  YV ) /DH 
AA22-( ( MASS-YVDOT ) * ( -MASS  *XG+NR ) - 
6  ( MASS*XG-NVDOT ) * ( -MASS+YR ) )/DH 

BBl 1- ( ( I Z -NRDOT ) * YDRS- ( MASS*  XG- YRDOT ) *NDRS )/DH 
BBl 2- ( ( IZ-NRDOT) * YDRB- (MASS *XG- YRDOT )*NDRB)/DH 
BB  2 1 - ( ( MASS - YVDOT ) *NDRS- ( MAS  S  *  XG-NVDOT ) *  YDRS ) /DH 
BB22- ( ( MASS-YVDOT ) *NDRB- ( MASS  * XG-NVDOT ) * YDRB )/DH 

WRITE  (*,1004) 
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READ  ( * , * ) 


RATIO 


BBl-BBll+RATIO*BBl2 
BB2-BB2 1-f  RATIO*  BB22 

EPS  -l.D-5 
I LHAX-1 500 

DO  1  I-1,ITC 

WRITE  (*,2001)  I , ITC 

TC-TCMIN+) 1-1 ) * ( TCMAX-TCMIN)/) ITC-1 ) 

TCD-TC*L/U 
POLEC-1 . 0/TCD 
ADI-3. 0*POLEC 
AD2  -  3 . 0  *  POLEC*  *  2 
AD3-POLEC*  *  3 
Al-BBl*U*U 
Bl-BB2*U*U 

Cl--ADl-(AAll+AA22)*U 
A2-) BBl *AA22-BB2*AAl 2  )  *U*  *  3 
B2  -  ( BB2  * AAl 1 -BB1 *AA2 1 ) *U*  *  3 
Kl-AD3/( ( BB2*AAll-BBl *AA21 ) *U*  *  3 ) 

C2-AD2-( AA11*AA22-AA12*AA21 ) *U* *2+BB2 *U*U*Kl 
K2-(Cl*B2-C2*Bl )/( Al *B2-A2*B1 ) 

K3-( C2  *A1-Cl *A2 )/( Al*B2-A2*Bl ) 

DO  2  J-1,IXD 

XD-XDMIN+ ( J-l ) * ( XDMAX-XDNIN )/( IXD-1 ) 

CALL  LINEAR) Kl , K2 , K3 , U , XD , AAl 1 , AAl 2 , AA21 , AA22 , BBl , BB2 , A, B, TL ) 
CALL  RGG(  5 , 5  ,  A,  B  ,  ALFR ,  ALF I  ,  BETA,  0  ,  ZZZ  ,  I  ERR  ) 

DO  11  IJE-1,5 

WR ( IJE)-ALFR( I JE )/BETA( IJE) 

WI( IJE)-ALFI ( I JE)/BETA( IJE) 

1  CONTINUE 

CALL  DSTABL( DEOS , WR,WI , FREQ) 

IF  (J.GT.l)  GO  TO  10 

DEOSOO-DEOS 

XDOO  -XD 

LL-0 

GO  TO  2 

0  DEOSNN-DEOS 

XDNN  -XD 
PR-DEOSNN*DEOSOO 
IF  (PR.GT.0.D0)  GO  TO  3 
LL-LL+1 

IF  (LL.GT.3)  STOP  1000 
I L-0 

XDO-XDOO 
XDN-XDNN 
DEOSO-DEOSOO 
DEOSN-DEOSNN 
6  XDL-XDO 

XDR-XDN 
DEOSL-DEOSO 
DEOSR-DEOSN 
XD- ( XDL+XDR ) /2 . DO 

CALL  LINEAR) Kl , K2 , K 3 , U , XD , AAl 1 , AAl 2 , AA21 , AA22 , BBl , BB2 , A, B , TL ) 
CALL  RGG) 5, 5, A, B, ALFR, ALF I , BETA , 0 , ZZZ , I  ERR ) 


DO  12  IJE-1,5 

WR ( I JE ) -ALFR ( I J E ) /BETA (  IJE) 
WX ( IJE) -ALFI ( IJE )/BETA(  IJE) 
12  CONTINUE 

CALL  DSTABL ( DEOS , WR , WI , FREQ ) 

DEOSM-DEOS 

XDM-XD 

P RL- D EOS L* DEOS M 

PRR-DEOSR*DEOSM 

IF  < PRL.GT. 0 . DO )  GO  TO  5 

XDO-XDL 

XDN-XDM 

DEOSO-DEOSL 

DEOSN-DEOSM 

IL-IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 
DIF-DABS(XDL-XDM) 

IF  (DIF. GT . EPS )  GO  TO  6 

XD-XDM 

GO  TO  4 

5  IF  ( PRR.GT. 0.D0)  STOP  3200 

XDO-XDM 
XDN-XDR 
DEOSO-DEOSM 
DEOSN-DEOSR 
IL-IL+1 

IF  ( IL.GT. ILMAX)  STOP  3100 
DI F-DABS ( XDM-XDR ) 

IF  (DIF. GT . EPS )  GO  TO  6 
XD-XDM 

4  LLL-10+LL 

WRITE  ( LLL , * )  XD/L.TC 
3  XDOO-XDNN 

DEOSOO-DEOSNN 
2  CONTINUE 
1  CONTINUE 


1001 

FORMAT 

(  ' 

ENTER 

MIN, 

MAX, 

AND 

INCREMENTS 

1002 

FORMAT 

(  ' 

ENTER 

MIN, 

MAX, 

AND 

INCREMENTS 

1003 

FORMAT 

<  ' 

ENTER 

U'  ) 

1004 

FORMAT 

(  ' 

ENTER 

BOW/STERN 

RUDDER  RATIO' ) 

1100 

FORMAT 

(  ' 

ENTER 

TIME 

LAG 

TL'  ) 

2001 

FORMAT 

(215) 

END 


C 


SUBROUTINE  DSTABL (DEOS, WR, WI , OMEGA) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  WR( 5) ,WI ( 5) 

DEOS--1 . 0D+20 
DO  1  1-1,5 

IT  (WR( I ) .LT.DEOS)  GO  TO  1 
DEOS-WR ( I ) 

I  J- 1 

1  CONTINUE 
OMEGA-WM  IJ) 

OMEGA-DABS ( OMEGA ) 

RETURN 

END 


C 


OF 

or 


TC'  ) 
XD'  ) 
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c 


SUBROUTINE  LI NEAR ( Kl , K2 , K 3 , U , XD , AA1 1 , AA12 ,  AA21 , AA2 2 , BBl , BB2 , 
&  A,  B , TL ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  Kl,K2,K3 
DIMENSION  A(5,5),B(5,5) 

A( 1 , 1 ) “0 . 0D0 
A(1,2)-0.0D0 
A( 1 , 3 ) “1 . 0D0 
A( 1 , 4 )-0 . 0D0 
A ( 1 , 5)-0.0D0 
A(  2 , 1 ) -0 . 0D0 
A( 2 , 2 } -0 . 0D0 
A( 2 , 3 )-0 . ODO 
A( 2 , 4 )-0 . 0D0 
A<2,5)-1 .ODO 

A( 3 , 1 )”( BB1*U*U*K1— BBl*U*U*U*Kl *TL/XD ) 

A( 3 , 2  j » ( AAl 1 *U+BBl *U*U*K2-BBl*U*U*Kl *TL/XD ) 

A ( 3 , 3 )-(AA12*U+BBl*U*U*K3+BBl *U*U*U* K1 *TL*TL/( 2 . D0*XD) ) 

A( 3, 4 )-(BBl*U*U*Kl/XD) 

A( 3 , 5  )  ■*  (  —  1 .DO+BBl*U*U*Kl *TL*TL/ ( 2 . DO  *  XD ) ) 

A(4,l)-U 
A( 4 , 2 )-l . ODO 
A( 4 , 3 ) -0 . ODO 
A( 4 , 4 )-0 . ODO 
A(  4 , 5 ) “0 . ODO 

A( 5, 1 )-(BB2*U*U*Kl-BB2*U*U*U*Kl*TL/XD) 

A( 5 , 2 )-( AA21 *U+BB2*U*U*K2-BB2*U*U*K1*TL/XD) 

A( 5 , 3 ) » ( AA22*U+BB2  *U*U*  K3  +  BB2  *U*U*U*K1 *TL*TL/( 2 . DO  *XD ) ) 
A(5,4)-(BB2*U*U*KJ/XD) 

A(5,5)-(BB2*  •;  'TL*TL/(  2  .  D0*XD )  ) 

B ( 1 , 1 ) -1 . ODO 
B ( 1 , 2 ) -0 . ODO 
B(l, 3)-0.0D0 
B( 1 , 4 )-0 . ODO 
B( 1 , 5 )“0 . ODO 
B ( 2 , 1 ) -0 . ODO 
B( 2 , 2 ) -I . ODO 
B( 2 , 3 ) -0 . ODO 
B ( 2 , 4 ) -0 . ODO 
B( 2 , 5 ) -0 . ODO 
B( 3 , 1 )-0 . ODO 
B ( 3 , 2 ) -0 . ODO 

B(  3, 3)-(BBl*U*U*U*Kl*TL*TL*TL/(6.D0*XD) ) 

B ( 3 , 4 ) »0 . ODO 

B(  3 , 5 ) -  (  BBl *U*U*Kl*TL*TL*TL/( 6 . DO*XD ) ) 

B ( 4 , 1 ) -0 . ODO 
B( 4, 2)-0.0D0 
B( 4, 3)-0.0D0 
B( 4 , 4 )-l . ODO 
B( 4 , 5 ) -0 . ODO 
B( 5 , 1 ) -0 . ODO 
B ( 5 , 2 )«0 . ODO 

B( 5, 3)-( 1 .0D0+BB2*U*U*U*Kl*TL*TL*TL/( 6.D0*XD) ) 

B( 5 , 4 ) -0 . ODO 

B( 5,5)-(BB2»U*U*Kl*TL*TL*TL/(6.D0*XD) ) 

RETURN 

END 
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APPENDIX  D 


C  PROGRAM  THESISN. FOR  (Tine  Delay  -  Newton's  Iteration  Method) 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  Kl , K2 , K3 , L , NR , NV, NDRS , NDRB , I Z , MASS , 

&  NRDOT , NVDOT 

DIMENSION  A( 4 , 4  )  , FVl ( 4 ) , IVl ( 4 ) ,ZZZ(4,4) ,  WR( 4 ) , WI ( 4 ) 

C 

OPEN  ( 11 , FILE- ' BIFl . RES ' , STATUS- ' NEW' ) 

C 

C  Vehicle  Parameters: 

WEIGHT-435.0 
IZ  -45.0 
L  -7.3 
RHO  -1.94 
G  -32.2 
XG  -0.0104 
MASS  -WEIGHT/G 
C 

YRDOT  — 0.00178*0. 5*RHO*L**4 
YVDOT  — 0.03430*0. 5*RHO*L**3 
YR  -+0.01187*0. 5*RHO*L**3 
YV  — 0.03896*0. 5*RHO*L**2 
YDRS  -+0. 02345*0. 5*RHO*L**2 
YDRB  -+0. 02345*0. 5*RHO*L**2 
NRDOT  — 0.00047*0. 5*RHO*L**5 
NVDOT  — 0.00178*0. 5*RHO*L**4 
NR  -  0.01022*0. 5*RHO*L**4 
NV  --0 . 00769*0 . 5*RHO*L**3 
NDRS  — 0. 337*0.02345*0.5*RHO*L**3 
NDRB  -+0. 283*0.02345*0.5*RHO*L**3 
C 

WRITE  (*,1001) 

READ  (*,*)  TCM I N , TCMAX , I TC 

WRITE  (*,1003) 

READ  ( *  ,  *  )  U 
C 

DH  -( IZ-NRDOT)*(MASS-YVDOT)- 
&  ( MASS*XG-YRDOT ) * ( MAS S*XG- NVDOT ) 

AAl l-( ( I Z -NRDOT) * YV- ( MASS*XG-YRDOT ) *NV)/DH 
AAl 2- ( ( I Z-NRDOT ) * ( -MASS+ YR ) - 
C  ( MASS  *XG-YRDOT) * ( -MASS*XG+NR ) )/DH 

AA2 1- ( ( MASS- YVDOT ) *NV- ( MASS*XG-NVDOT ) *  YV ) /DH 
AA22-( ( MASS- YVDOT) * < -MASS *XG+NR ) - 
&  ( MASS*XG-NVDOT ) * ( -MASS+YR ) )/DH 

BBl 1— ( ( I Z-NRDOT ) * YDRS- ( MASS*XG-YRDOT ) *NDRS )/DH 
BBl 2- ( ( I Z-NRDOT) *YDRB-( MASS*XG-YRDOT) *NDRB )/DH 
BB  2 1 - ( ( MASS- YVDOT ) *NDRS - ( MASS  *  XG-NVDOT ) *  YDRS ) /DH 
BB2  2- ( ( MASS- YVDOT ) *NDRB- ( MASS  *XG-NVDOT ) *  YDRB ) /DH 
C 

WRITE  (*,1004) 

READ  {*,*)  RATIO 

C 

BBl-BBl l+RATIO*BBl 2 
BB2-BB2 1  +  RATIO*  BB2  2 

C 

WRITE  (*,1005) 

READ  ( *, * )  TL 

C 

EPS-1 .D-10 
ITMAX-1000 
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non  non 


c 

DO  1  I-1,ITC 

WRITE  (*,2001)  I , ITC 

TC-TCMIN+ ( 1-1 ) * ( TCMAX-TCMIN )/( ITC-1) 

TCD-TC*L/U 
POLEC-1 .0/TCD 
ADl-3 . 0*POLEC 
AD2- 3 . 0*POLEC*  *2 
AD3-P0LEC*  *  3 
Al-BB1*U*U 

Bl-BB2*U*U  * 

Cl— ADl-(AAll+AA22) *U 
A2-(  BB1*AA22-BB2*AA12 )  *11**3 
B2- ( BB2 *AAl 1-BBl *AA21 ) *U**3 
Kl«AD3/( (BB2*AA11-BB1*AA21)*U**3) 

C2-AD2-( AAll*AA22-AA12*AA21 ) *U* * 2+BB2 *U*U*Kl 
K2-(C1*B2-C2*B1 )/( A1 *B2-A2*B1 ) 

K3-(  C2*A1-C1 *A2 )/( Al*B2-A2  *B1 ) 

Al-  ( AAl 1 *BB2-AA21 *BBl )*U*U*U*Kl 

A2-  (AAl 1 *AA22-AAl2 *AA21 ) *U*U+ ( AAl 1  * BB2-AA21 *BBl )*U*U*U*K3+ 

&  ( AA22  *BBl -AAl 2  *BB2 ) *U*U*U*K2-BB2*U*U*Kl 

A3— ( AAl 1+AA22 ) *U-( BBl*K2+BB2*K3 )  *U*U 
BO-  (AAll*BB2-AA21*BBl)*U*U*U*U*Kl 
B1-- ( AAl 2  *BB2-AA22  *BBl ) *U*U*U*Kl-BB2*U*U*U*Rl 
B2  —  BBl*U*U*Kl 

COMPUTE  INITIAL  APPROXIMATION  FOR  OMEGA  (TL-0) 

AQ-B2*A3-B1 
BQ-Bl *A2-B2  *Al-B0*A3 
CQ-B0*A1 

DET-BQ*  *2-4 . DO  * AO *CQ 
IF  (DET.LT.0.D0)  GO  TO  1 
ZTl- ( -BQ+DSQRT ( DET ) )/(2.0D0*AQ) 

ZT2- ( -BQ-DSQRT( DET ) )/(2.0D0*AQ) 

IF  (ZTl .GE.0.D0)  OM-DSQRT(ZTl) 

IF  (ZT2.GE.0.D0)  OM-DSQRT( ZT2 ) 

ALPHAl-AQ 
ALPHA2-BQ 
ALPHA3-CQ 
BETAl  — B2 

BETA2-B0+B2*A2-Bl*A3 
BETA3-Bl*Al-B0*A2 

COMPUTE  EXACT  OMEGA  USING  NEWTON'S  ITERATION 
OMOLD-OM 

3  F- ( BETAl *OMOLD* *  5  +  BETA2 *OMOLD*  * 3+BETA3 *OMOLD) *DSIN(OMOLD*TL ) 

&  +(ALPHAi*OMOLD** 4 +ALPHA2*OMOLD* *  2  + ALPHAS ) *DCOS ( OMOLD*TL ) 

F PR I ME- ( 5 . D0*BETA1 *OMOLD*  *  4  +  3 . DO  *BETA2  *OMOLD*  *  2  + BETA 3 ) 

&  *DS IN ( OMOLD*TL )  +  ( BETAl *OMOLD*  *  5  +  BETA2 *OMOLD*  *  3  +  BETA3 *OHOLD ) 

&  *TL*DCOS(OMOLD*TL )  +  ( 4 . D0*ALPHAl *OMOLD* *  3  +  2 . D0*ALPHA2  *OMOLD )  * 

t  *  DCOS <  OMOLD*TL ) - ( ALPHAl *OMOLD*  *  4+ALPHA2  *OMOLD* *  2+ALPHA3 ) 

&  *TL*DSIN( OMOLD*TL ) 

IF  ( FPRIME . EQ . 0 . DO )  STOP  1112 
IF  (F.EQ.O.DO)  GO  TO  2 
OMNEW-OMOLD-F/FPRIME 
OMD I F-DABS ( OMNEW-OMOLD ) 

IF  (OMDIF.LE.EPS)  GO  TO  2 
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OMOLD-OMNEW 

IT-IT+1 

IF  ( IT.GT. ITMAX)  STOP  1111 
GO  TO  3 
C 

2  OM-OMNEW 

XDNUM-DSQRT( ( BO-B2*OM* *  2 ) *  *  2+ ( Bl *OM ) * *2 ) 

XDDEN-OM*  DSQRT( ( A1-A3  *OM*  *  2 ) *  * 2+ ( A2  *OM-OM*  *  3 ) *  *  2 ) 
XD-XDNUM/XDDEN 
XD-XD/L 
C 

WRITE (11,*)  XD , TC , OM 
C 

1  CONTINUE 

1001  FORMAT  ('  ENTER  MIN,  MAX,  AND  INCREMENTS  OF  TC' ) 

1003  FORMAT  ('  ENTER  U') 

1004  FORMAT  ('  ENTER  BOW/STERN  RUDDER  RATIO') 

1005  FORMAT  ('  ENTER  TIME  LAG') 

2001  FORMAT  (215) 

END 
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